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Summary. We report numerical simulations of strongly vibrated granular materials 
designed to mimic recent experiments performed both in presence [1] or absence 
[2] of gravity. We show that a model with impact velocity dependent restitution 
coefficient is necessary to bring the simulations into agreement with experiments. 
We measure the scaling exponents of the granular temperature, collision frequency, 
impulse and pressure with the vibrating piston velocity. As the system changes 
from a homogeneous gas state at low density to a clustered state at high density, 
these exponents are all found to decrease continuously with the particle number. In 
absence of gravity, a loose cluster appears near the upper wall, opposite the piston, 
and acts as a buffer for fastest particles leading to unexpected non-extensive scaling 
exponents ; whereas in presence of gravity, the cluster bounces as a single inelastic 
body. All these results differ significantly from classical inelastic hard sphere kinetic 
theory and previous simulations, both based on a constant restitution coefficient. 



1 Introduction 

Granular gases [3] often seem to be the exclusive domain of theory and sim- 
ulation. For example, the "freely coohng" granular gas [3] is impossible to 
experimentally realize: even if an experimentalist somehow managed to im- 
plement periodic boundaries in all directions, he would still have to contend 
with energy and time scales that vary over many orders of magnitude. Another 
frequently studied example of granular gases is the vibrated granular medium. 
This case is much closer to experiments, because it is relatively simple to put 
particles in a box and shake them, or drive particles with a vibrating piston. 
Indeed, there are many experimental [4-6], numerical [6-8], and theoretical 
[4,9-11] studies of this system. However, numerous questions remain about 
the link between experiments on one hand, and theory and simulations on 
the other. Most numerical and theoretical studies were not intended to be 
compared with experiments. Therefore, they have parameter values far from 
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the experimental ones, and none of them predict even the most basic features 
of the experimental results. For example, numerical and theoretical studies 
often assume that the vibration frequency is very high so that the shaking 
wall can be replaced by a thermal boundary. The vibrated granular gas thus 
attains a steady state. But very high frequencies with an amplitude sufficient 
to fluidize the granular materials are difScult to attain experimentally; the 
vibration amplitude in experiments is often a significant fraction of the size 
of the container. In experiments, the pressure and granular temperature all 
vary strongly over the period of one vibration cycle. It has therefore been im- 
possible to compare the experimental and numerical results in a meaningful 
way. 

In this paper, we bridge the gap between experiments and numerics by pre- 
senting simulations of strongly vibrated granular materials designed to mimic 
recent experiments [1,2]. We present the first simulations which resemble the 
experiments for a large range of parameters. We show that two parameters are 
especially important for the agreement between experiment and simulation. 
First of all; an impact velocity dependent restitution coefficient is necessary to 
bring the simulations into agreement with experiment. Most previous studies 
consider only constant restitution coefficient [6-8]. Secondly, it is important to 
explicitly consider the number of particles N. Studying only one value of TV or 
comparing results obtained at different N can lead to interpretive difficulties. 

We measure the scaling exponents of the granular temperature, collision 
frequency, impulse and pressure with the vibrating piston velocity, and find 
results that differ significantly from classical inelastic hard sphere kinetic the- 
ory and previous simulations. Finally, we show that the system undergoes a 
smooth transition from a homogeneous gas state at low density to a clustered 
state at high density. Depending on whether gravity is present or not, the tran- 
sition manifests itself in different ways: due to inelastic collisions (through the 
restitution coefficient properties) , a nearly motionless loose cluster is observed 
in microgravity with non-extensive scaling properties, whereas a bouncing 
cluster appears under gravity. 

2 Description of the simulations 
2.1 The Vciriable coefficient of restitution 

The greatest difference between our simulations and the previous numerical 
studies of vibrated granular media [6-8] is that we use a restitution coefficient 
that depends on impact velocity. The restitution coefficient r is the ratio 
between the relative normal velocities before and after impact. In all previous 
simulations of strongly vibrated granular media, the coefficient of restitution 
is considered to be constant and lower than 1. However, since a century, it 
has been shown from impact experiments that r is a function of the impact 
velocity v [12-16]. Indeed, for metallic particles, when v is large {v > 5 m/s 
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[13]), the colliding particles deform fully plastically and r oc v~^/^ [12-14]. 
When w < 0.1 m/s [13], the deformations arc clastic with mainly viscoclastic 
dissipation, and 1 — r oc v^/^ [14-17]. Such velocity-dependent restitution 
coefficient models have recently shown to be important in numerical [18- 
23] and experimental [16, 24] studies. Applications include; granular fluidlike 
properties (convection [19], surface waves [20]), collective coUisional processes 
(energy transmission [21], absence of collapse [16,22]) and planetary rings 
[23, 24]. But surprisingly, such model has not yet been tested numerically for 
strongly vibrated granular media. 
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Fig. 1. The restitution coefficient r as a function of impact velocity v, as given in 
Eq. (1) (solid line). The dashed lines show «o = 0.3 m/s and ro = 0.95. Experimental 
points (•) for steel spheres were extracted from Fig.l of Ref [25] 



In this paper, we use a velocity dependent restitution coefficient r{v) and 
join the two regimes of dissipation (viscoelastic and plastic) together as simply 
as possible, assuming that 



r{v) = 



1 
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where Wg = 0.3 m/s is chosen, throughout the paper, to be the yielding velocity 
for stainless steel particles [13,25] for which rg is close to 0.95 [25]. Note that 
'^0 ~ where p is the density of the sphere [13]. We display in Fig. 1 

the velocity dependent restitution coefficient of Eq. (1), with ro = 0.95 and 
vo = 0.3 m/s, that agrees well with experimental results on steel spheres from 
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Ref. [25]. As also already noted by Ref. [13], the impact velocity to cause 
yield in metal surfaces is indeed relatively small. For metal, it mainly comes 
from the low yield stress value {Y ^ 10^ N/m^) with repect to the elastic 
the Young modulus {E ~ 10^^ N/m^). Most impacts between metallic bodies 
thus involve some plastic deformation. 

2.2 The other simulation parameters 

The numerical simulation consists of an ensemble of identical hard disks of 
mass m « 3 X 10~^ kg excited vertically by a piston in a two-dimensional box. 
Simulations are done both in the presence and absence of uniform gravity g. 
Collisions are assumed instantaneous and thus only binary collisions occur. 
For simplicity, we neglect the rotational degree of freedom. Collisions with the 
wall are treated in the same way as collisions between particles, except the 
wall has infinite mass. 

Motivated by recent 3D experiments on stainless steel spheres, 2 mm in 
diameter, fluidizcd by a vibrating piston [1], we choose the simulation pa- 
rameters to match the experimental ones: in the simulations, the vibrated 
piston at the bottom of the box has amplitude A = 25 mm and frequencies 
5 Hz < / < 50 Hz. The piston is nearly sinusoidally vibrated with a wave- 
form made by joining two parabolas together, leading to a maximum piston 
velocity given by V = 'iAf. The particles are disks d = 2 mm in diameter 
with stainless steel collision properties through vq and ro (see Fig. 1). The 
box has width L = 20 cm and horizontal periodic boundary conditions. Since 
our simulations are two dimensional, we consider the simulation geometri- 
cally equivalent to the experiment when their number of layers of particles, 
n = Nd/L, are equal. Hence in the simulation, a layer of particles, n = 1, 
corresponds to 100 particles. We checked that n is an appropriate way to 
measure the number of particles by also running simulations at L = 10 cm 
and L = 40 cm. None of this paper's results depend significantly on L, excpet 
in Sec. 4.2, where this point is discussed. As in the experiments, the height 
h of the box depends on the number of particles in order to have a constant 
difference h — ho = 15 mm, where ho is the height of the bed of particles at 
rest. Heights are defined from the piston at its highest position. 

3 Comparison of simulation and experiment 

3.1 The importance of the Vciriable coefficent of restitution 

We examine first the dependence of the pressure on the number of particle 
layers for maximum velocity of the piston l<y<5m/s(F = 4^4/). The 
time averaged pressure at the upper wall is displayed in Fig. 2 as a function of 
n for various /: from the cxpc^riments of Fak;on et al. [1] (sec Fig. 2a), from our 
simulations with velocity dependent restitution coefficient r = r(v) proposed 
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in Eq. (1) (see Fig. 2b), and with constant restitution coefficient r = 0.95, often 
used to describe steel particles (see Fig. 2c), and finally with an unrealistic 
constant restitution coefficient r = 0.7 (see Fig. 2d). Simulations with r = 
r{v) give results in agreement with the experiments: At constant external 
driving, the pressure in both Figs. 2a and 2b passes through a maximum for a 
critical value of n roughly corresponding to one particle layer. For n < 1, most 
particles are in vertical ballistic motion between the piston and the lid. Thus, 
the mean pressure increases roughly proportionally to n. When n is increased 
such that n > 1, interparticle collisions become more frequent. The energy 
dissipation is increased and thus the pressure decreases. As we show below in 
Sec. 4.1, this maximum pressure is not due to gravity because it also appears 
in simulations with g = and r = r{v). Turning our attention to Fig. 2c, 
we see that setting r = 0.95 independently of impact velocity gives pressure 
qualitatively different from experiments. The difference between Figs. 2b and 
2c can be understood by considering a high velocity collision (e.g. v = 1 m/s). 
In Fig. 2b, this collision has a restitution coefficient of r = r(l m/s) ~ 0.7 (see 
Fig. 1), whereas in Fig. 2c, r is fixed at 0.95 for all collisions. This means that 
for equal collision frequencies, dissipation is much stronger for r = r{v) than 
for r = 0.95, because the high velocity collisions dominate the dissipation. 
Stronger dissipation leads to lower granular temperatures and thus to lower 
pressures. 

We can check this interpretation by changing the constant restitution co- 
efficient to r = 0.7 and then comparing it to r = r(f). In these two cases, 
the high velocity collisions will have roughly the same restitution coefficient. 
We indeed observed a pressure that decreases for large n for constant r = 0.7 
(see Fig. 2d). Therefore, surprisingly, constant r = 0.7 reproduces more pre- 
cisely the experimental pressure measurements than constant r = 0.95, even 
though r = 0.95 or r = 0.9 is often given as the restitution coefficient of 
steel. Consulting Fig. 1, one can sec that r{v) = 0.7 is well inside the plastic 
regime. Therefore, our results suggest that plasticity is much more important 
than visco-elasticity in vibrated granular media. Although constant r = 0.7 
gives good agreement with the experiments in Fig. 2, this will not be true in 
all situations. If the inter-particle collisions were much gentler, the effective 
coefficient of restitution would rise, and another choice of constant r would 
be needed. Furthermore, the pressure is just one property of the system. If we 
look at other property, we see that r = 0.7 does not give good results. An ex- 
ample is shown in Fig. 3, where we show three snapshots from three different 
simulations with n = 3 in absence of gravity. When r = r{v), the particles 
are concentrated in the upper half of the chamber, but they are evenly spread 
in the horizontal direction (see Fig. 3a). The system is hotter and less dense 
near the vibrating wall, and colder and denser by the opposite wall. Such a 
loose cluster has been already observed in microgravity during parabolic flight 
experiments [26], where grains are agitated with a piston. In contrast, when 
r = 0.95, the particles arc imiformly distributed throughout the domain (see 
Fig. 3b). This occurs because the particle move much faster than the piston. 
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(a) experiments (b) r = r{v) 




(c) r = 0.95 (d) r = 0.7 
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Fig. 2. Time averaged pressure P on the top of tlie cell as a function of particle 
layer, n, for various vibration frequencies, / : (a) Experimental results from [1] for 
stainless steel beads 2 mm in diameter, with ^4 = 25 mm, 10 < / < 2QHz with a 1 
Hz step (from lower to upper) and h — ho — 5 mm. (b) Numerical simulation where 
the coefficient of restitution is given by Eq. (1), (c) Numerical simulation with 
a coefficient of restitution is 0.95, independent of impact velocity, (d) Numerical 
simulation with a coefficient of restitution is 0.7, independent of impact velocity. 
The simulations (b) - (d) are 2D with gravity, done for 2 mm disks, with A = 25 
mm, 10 < / < 30 Hz with a 2 Hz step (from lower to upper) and h — ho = 15 mm. In 
the simulations, the two-dimensional pressure is given in units of mg/d « 0.15 Pa-m. 

and thus can fill the space left by the descending piston. Finally, when r = 0.7, 
the majority of the particles are confined to a tight cluster, pressed against the 
upper wall, coexisting with low density region (see Fig. 3c). This instability 
has been already been reported numerically [27], although for much differ- 
ent parameters (constant restitution coefficient r = 0.96, thermal walls, no 
gravity, and large n). However, nothing like this was ever seen experimentally. 
Therefore, if one is seeking information about particle positions, r = 0.7 gives 




Fig. 3. Snapshots from the simulations with n = 3, gravity p 7^ 0, driving frequency 
/ = 30 Hz and h — ho = 15 mm. The upper wall is stationary, and the lower wall 
is the piston, which is at its lowest position in all three snapshots. The horizontal 
boundaries arc periodic (indicated by dashed lines). Gravity points downwards, (a) 
r = r{v), as given in Eq. (1), (b) constant r = 0.95, and (c) constant r = 0.7. In 
(c) we see a tight cluster which was not observed in the experiments. 



incorrect results even though it gives acceptable results for the pressure. We 
conclude, therefore, that the only way to successfully describe all the prop- 
erties in all situations is to use a velocity dependent restitution coefficient 
model. 



3.2 The importance of the particle number 

Observations suggest that the the pressure (or granular temperature) of 
strongly vibrated granular medium is controlled by the piston vibration veloc- 
ity V. It is important to point out that V is not the only way to characterize 
the vibration. One could also use the maximum piston acceleration 7^. When 
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Fig. 4. The exponents 8 as a function of n wiiicii give tlic scaling of the granular 
temperature T (0), collision frequency Cupper (*), mean impulsions AI (A) and 
pressure P (o). All these quantities are proportional to V^^"K Without gravity; (a) 
for r = 0.95 and (b) for r = r{v). With gravity: (c) for r = 0.95 and (d) for r = r{v). 
The exponents are obtained by fixing n and performing eleven simulations, varying 
/ from 10 Hz to 30 Hz. Then logX (where X is the quantity being considered) is 
plotted against log V. The resulting curve is always nearly a straight line (except for 
n > 3 in (d) - see text), and the exponent is calculated from a least squares fit. The 
pressure scaling point (•) on (b) is from experiment [2] performed in low-gravity. 
See Fig. 3a (resp. Fig. 3b) for typical snapshots corresponding to n = 3, p and 
r = r{v) (resp. r = 0.95) 



r is close to g, it controls the behavior of the system, i.e., adjusting A and / 
while keeping F constant does not change the system's behavior much. But in 
the simulations presented here, F ^ g, and the system's behavior is controlled 
by V. This can be checked by multiplying the frequency by 10 while dividing 
A by 10, thus keeping V the same (while F increases by an order of mag- 
nitude). Doing so changes the pressure only by about 20%. It is reasonable, 
therefore, to look for scaling relations in V. 

Many authors have postulated that the pressure on the upper wall P (or 
granular temperature T) is related to the piston velocity V through P,T oc 
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V^. However, it is not clear what the correct "scaling exponent" 6 should 
be. This question has been addressed several times in the past, without a 
clear resolution of the question [8-11]. For example, kinetic theory [4, 10] and 
hydrodynamic models [11] predict T oc V'^ whereas numerical simulations 
[6, 7] or experiments [4-6] give T oc , with 1 < 9 < 2. These studies were 
done at single values of n. In this section, we show that it is very important 
to explicitly consider the dependence of the scaling exponents on n. We also 
consider the effect of gravity and a variable coefficient of restitution. Doing so 
enables us to explain and unify all previous works. 

At the upper wall, we measured numerically the collision frequency. Cupper 
and the mean impulsion per collision, AI for various frequencies of the vibrat- 
ing wall and numbers of particles in the box, with r = r{v) or with r = 0.95, 
in the presence or absence of gravity. The time averaged pressure on the upper 
wall can be calculated from these quantities using 

P = CupperAI/L . (2) 

(By conservation of momentum, the time averaged pressure on the lower wall 
is just P plus the weight of the particles Nmg.) The total kinetic energy of 
the system is also measured to have access to the granular temperature, T. 
Cupper, AI^ P, and T are all found to fit with power laws in for our range 
of piston velocities. Fig. 4 shows exponents of Cupper, AI, P, and T as a 
function of n. When g = and r is constant (see Fig. 4a), we have P ^ V"^, 
AI ^ V and Cupper ~ V for all n. We call these relations the classical ki- 
netic theory scaling. This scaling can be established by simple dimensional 
analysis when the vibration velocity V provides the only time scale in the 
system. This is the case for g = and r independent of velocity. However, in 
the experiments, two additional time scales are provided, one by gravity and 
another by velocity dependent restitution coefficient. Numerical simulations 
can separate the effects of these two new time scales on the scaling exponents 
6. This is done in Fig. 4b [where ,9 = but r = r{v)] and Fig. 4c [where r is 
constant but (? ^ 0]. In both figures, all the exponents become functions of 
n. However, the time scale linked to r = r{v) leads to much more dramatic 
departure from the classical scaling. After considering the two time scale sepa- 
rated, let us consider the case corresponding to most experiments, where both 
gravity and r — r{v) are present (see Fig. 4d). The similarity between this 
figure and Fig. 4b confirms that the velocity dependent restitution coefficient 
has a more important effect than the gravity. Furthermore, only the variation 
of the restitution coefficient on the particle velocity explains the experiment 
performed in low-gravity [2]. This experiment gives a V^^^ pressure scaling 
(•-mark on Fig. 4b) for n = 1 and a motionless clustered state for n > 2. 
Only the simulation with r = r{'v) can reproduce these results (see Fig. 4b) 
whereas constant r simulations leads to the classical scaling {P ocV^, see Fig. 
4a) and only a gaseous state for all n shown in the figure. 

As shown in Fig 4, it is thus very important to explicitly consider the 
dependence of 9 on n. In all cases, except the unrealistic case of Fig. 4a, 9 
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depends on n. To our knowledge, the only experiment [1] to systematically 
investigate this effect shows that T oc V^^"\ with 6 continuously varying 
from 6 = 2 when n ^ 0, as expected from kinetic theory, to ~ for large 
n due to the clustering instability. These experiments [1] performed under 
gravity (shown in Fig. 5) arc well reproduced by the simulations of Fig. 4d. In 
both cases, the observed pressure and granular temperature scaling exponents 
strongly decrease with n. 
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Fig. 5. Experimental data performed under gravity from Ref. [1]: The exponents 
6{n) of time averaged pressure (□) (see Fig. 2a), and kinetic energy extracted from 
density profile (o) or volume expansion (0) measurements. These data should be 
compared with the simulations of Fig. 4d. 



We finish this section by noting two curious facts about Fig. 4. First of 
all, in Fig. 4b [5 = and r = r{v)], 6 « 1 for the pressure and temperature 

when n > 2. This is the sign of new robust scaling regime where P, T oc V, 
which is the subject of Sec. 4.1. Secondly, in Fig. 4d [g 7^ and r = r{v)], 
the scaling exponents are not shown for n > 3, because the dependence of P, 
T Copper and AI on V is no longer a simple power law. (More precisely, we 
do not plot a point on Fig. 4 when |logXobserved — logXfittedl > 0.25 for any 
of the eleven simulations used to calculate the exponent - see caption.) The 
power law breaks down because there is a resonance between the time of flight 
of the cluster under gravity and the vibration period. This is the subject of 
Sec. 4.2. 
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4 Effects of clustering (n > 2 or 3) 
4.1 Clustering without gravity 

In this section, we examine the situation concerning reaUstic particles [r = 
r{v)] in microgravity (g = 0). At high enough density, a loose cluster is formed 
that resembles Fig. 3a (even though g =/= in Fig. 3a), and we noted that 
P (X V for n > 2 (see Fig. 4b). We would like to know if there is some simple 
physical explanation for this aparently robust scaling relation. To investigate 
this question, we present in Fig. 6 the dependence of PonV,n and h — ho- 



(a) h- ho = 15 mm (b) / = 20 Hz 




n 



Fig. 6. Time averaged pressure P on the top of the cell as a function of particle 

layer, n: (a) for various vibration frequencies: / = 10 Hz (lowest curve) to / — 30 Hz 
(uppermost curve) with steps of 2 Hz; (b) for various heights: h — ho — 5 mm 
(uppermost curve) to h — ho = 50 mm (lowest curve) with steps of 5 mm. Inset: P 
vs. {h — ho)~^ for n = 10 and 10 cm < h — ho < 50 cm. Simulations were done with 
g = and r = r{v). 

Fig. 6a shows the same qiiantities as Fig. 2, except that the scale of the .r-axis 
has changed: instead of < n < 4.5, we show < n < 10. Note that Fig. 6a 
concerns the same situation as in Fig. 2b, except here the gravity has been 
"turned off". Thus it is not surprising that for n < 3, the two figures are 
almost the same: after a rapid increase of the pressure for small n, there is 
a maximum pressure near n w 1. But for n > 3, the figures are different: 
in Fig. 2b, the pressure decreases as particles are added, but in Fig. 6a the 
pressure is nearly independent of n. Therefore, when g = 0, the pressure is 
not changed by adding more particles, as long as there are already enough 
particles present. Considering now the dependence of P on V, we note that 
the evenly spaced curves in Fig. 6a suggest that P oc V. We have confirmed 
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this proportionality by calculating the pressure scaling exponent and checking 
that it remains close to 1. (The exponent for n < 6 has already been shown 
in Fig. 4b.) In Fig. 6b, we examine how changing h — ho affects the pressure. 
We see that the pressure is independent of n for n > 3 and h — ho> 10, with 
the pressure decreasing as the height increases. Examining the dependence of 
the pressure onh — ho shows that P (x (h — ho)~^. (For example, see the inset 
in Fig. 6b.) Thus, the pressure obeys the simple non-extensive relation 



Pa 



h — ho 



(3) 



Recall that this relation concerns high enough number of realistic particles 
in microgravity. Therefore, this pressure scaling of Eq. (3) may be observable 
in microgravity experiments. One complication is that in microgravity exper- 
iments, one does not usually agitate the grains with a piston (except very 
recently [26]) - it is much simpler to shake a box of fixed size [2]. 



(a) h — ho = 15 mm 



(b) / = 20 Hz 




Fig. 7. The particle-piston collision rate Ciower as a function of particle layer, n: (a) 
for various vibration frequencies: / = 10 Hz (lowest curve) to / = 30 Hz (uppermost 
curve) with steps of 2 Hz; (b) for various heights: h—ho = 10 mm (uppermost curve) 
to h — ho = 50 mm (lowest curve) with steps of 5 mm. Simulations are done with 
g = and r = r{v). 



In Fig. 3a, we observe that the majority of particles remain in a loose 
cluster pushed against the stationary wall, opposite the piston. Only those 
particles that "evaporate" from the cluster are struck by the piston. The 
evaporation rate can be estimated from the rate of collisions between the 
piston and the particles Ciowor- This collision rate has a very curious behavior, 
as shown in Fig. 7. As can be seen from Fig. 7a, Ciower is roughly independent 
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of n and V when n > 3. This behavior hold for other values of h — ho, as can 
be seen from Fig. 7b. Thus, at high enough density, the collision frequency on 
the vibrating wall is found to be 

Clower OC -— . (4) 

II — llQ 

We now present a simple model and a physical interpretation of the unex- 
pected scaling of Eq. (3) for such a loose cluster. We suppose that the cluster 
against the stationary wall has a granular temperature independent of both 
the number of particles N and the vibration velocity V. The granular tem- 
perature is independent of N and V because the impact velocity dependent 
restitution coefficient acts as a kind of "granular thermostat" . As one can see 
from Eq. (1) or Fig. 1, when the typical particle velocity v is greater than 
vo, r{v) is small, and energy is dissipated rapidly. On the other hand, when 
V < vq, r{v) is nearly unity and energy dissipation is slow. Thus, in a weakly 
forced cluster, like the one found at the stationary wall, T ~ will be near 
Vq, independent of the strength of the forcing and the size of the cluster. 

Since the granular temperature of the cluster is independent of N and 
V, so is the flux Fevap of grains "evaporating" from the cluster and moving 
towards the piston. If all the grains that evaporate from the cluster reach 
the piston, then we have Ciowcr = ^ovap N'^V'^. In their collision with the 
piston, the grains acquire an upwards velocity proportional to V. Thus the 
flux of momentum entering the system at the piston is FCiower oc N^V^. Since 
momentum is conserved, the flux of momentum leaving the system through 
the stationary wall must have the same value. But the time averaged pressure 
P is just the time averaged momentum flux divided by the area of the upper 
wall. Thus our model predicts P oc iV^V^-'^, in agreement with Eq. (3). 

This theory also explains why P and Ciower are inversely proportional to 
h — ho- When a particle evaporates from the cluster, it must travel a certain 
distance before it encouters the piston. This distance increases with h — ha and 
thus the particle's travel time also increases. During its voyage, the evaporated 
particle could be struck by another particle that has just encountered the 
piston. If this happens, both particles are scattered back into the cluster. 
Thus the evaporated particle never reaches the piston. If we assume that the 
probability of an evaporated particle being scattered back into the cluster is 
independent of time, the number of particles reaching the piston is inversely 
proportional to h — ho. 

This explanation can be checked by looking at the probability density of 
vertical {y) velocities parallel to the direction of vibration. In Fig. 8, we show 
these density functions. Note the assymmetry of the distibution's positive 
and negative wings. The edge of the positive wing increases linearly with V 
whereas the negative wing extends much more slowly. The reason for this is 
that the positive wing is populated by particles which have just been struck 
by the piston. Their velocity is thus proportional to V. On the other hand, 
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Fig. 8. The probability density function of the vertical velocities at high density, 
for n = 5, r = r{v), and 5 = 0. The driving frequencies are / = 10 Hz (narrowest 
distribution) to / = 50 Hz (broadest distribution) with steps of 10 Hz. 

the negative wing is due to fluetuations in the chister. The negative wing thus 
expands slowly because of the "thermostat" effect described above. 

We conclude by emphasizing that this effect can only be seen with the 
velocity dependent restitution coefficient. If one uses a constant restitution 
coefficient, one obtains P oc y^, no matter the value of r. 

4.2 Clustering with gravity 

In this section, we consider the situation that applies to most of the experi- 
ments, i.e. realistic particles [r = r{v)] under gravity. We will show that the 
relatively simple situation described in the previous section is complicated 
by the presence of gravity. Under gravity, the loose cluster will not remain 
near the upper wall (unless / is large), but will fall downwards towards the 
piston. At certain values of /, h — ho and large enough n, a resonance occurs 
between the driving frequency and the time of flight of the particles, leading 
to a cluster that bounces back and forth between the piston and the wall. As 
/ increases, however, we recover the situation found in the previous section. 

In Fig. 9, we show two pictures from the experiments of Ref. [1]. At low 
densities, the particles move like atoms in a gas (Fig. 9a). At higher densities, a 
cluster forms that bounces up and down on the piston (Fig. 9b). This cluster 
was observed only for certain frequencies, densities and container heights. 
It was not possible to establish a criterion for the existence of the cluster 
based on the experimental data. Both of these behaviors are reproduced in 
the simulations (Fig. 10). Therefore, the simulations should give insight into 
the origin of the bouncing cluster. 





Fig. 9. Photographs of particle trajectories from the experiments of [1], showing the 
cluster at high densities. In both pictures f — 20 Hz and A = 40 mm. The number 
of layers of particles is (a) n = 1 (b) n = 4. The driving piston is at the bottom 
(not visible), the inner diameter of the tube being 52 mm. 

(a) n = 1 (b) n = 5 




Fig. 10. Particle trajectories from the simulations with / = 20 Hz, A = 25 mm, 

r = r{v), and h — ho — 15 mm. The number of layers of particles is (a) n = 1 and 
(b) n = 5. The containers have different sizes, because the size of the container is 
increased as more particles are added, as was done in the experiments. 
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We first examine the relation between P and F in a case where the bounc- 
ing cluster is observed. One docs not obtain the simple scaling in Eq. (3) for 
the pressure. Instead, one obtains a much more complicated relation between 
P and V. Indeed, as we show in Fig. 11, the pressure displays two maxima 
as a function of V. These maxima correspond to the bouncing cluster shown 




V (m/sec) 



Fig. 11. Pressure as a function of the piston velocity, V = 4Af, for n = 6, r = r{v), 
h — ho = 15 mm under gravity. Simulations where done for / = 5 Hz to / = 50 Hz 
with steps of 0.2 Hz. The dots mark the particular frequencies shown in Fig. 12. 

in Figs 9b and 10b. The bouncing cluster transports momentum very effec- 
tively between the piston and the wall. When the bouncing cluster disappears 
at n < 3, the maxima also disappear. These maxima make it impossible to 
describe the relation between the pressure and the forcing by a simple power 
law. This is the reason why there are no points on Fig. 4d for n > 3. 

Unlike the other results of this paper, Fig. 11 changes significantly when 
L and N are changed while holding n fixed. When L and are both divided 
by two, the first maximum moves to F w 1 m/sec and the second maxi- 
mum disappears. When L and N are doubled, the second maximum moves 
to near V « 1.4 m/sec. Therefore, the results of this section should not be 
considered as a systematic investigation of the resonance, but rather a de- 
scriptive introduction. The resonance is sensitive to L because it can involve 
horizontal collective motions of the particles, not just vertical ones. Indeed, 
for L = 40 cm, one can observe waves similar to those discussed in [20] . 

To investigate the origins of the bouncing cluster, we show in Fig. 12 
the position of the piston and the approximate location of the cluster as a 
function of time for various driving frequencies. The frequencies are chosen to 
correspond to the maxima and minima of the pressure, and arc indicated by 
the large points in Fig. 11. At very low frequencies (Fig. 12a), the pressure 
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(a) / = 5 Hz, V = 0.5 m/sec (b) / = 8.2 Hz, V = 0.82 m/sec 




(c) / = 13.8 Hz, y = 1.38 m/scc (d) / = 21.6 Hz, v = 2.16 m/scc 




(e) / = 26.4 Hz, V = 2.64 m/soc (f) / = 50 Hz, "1/ = 5 m/sec 

5r ' 1 ' 1 r ' 1 ' 



4- 




D.O 0.1 0.2 0.1 0.2 



time (sec) time (sec) 

Fig. 12. The motion of the grains under gravity at various vibration frequencies 
during 0.2 seconds, ior n = 6, h — ho = 25 mm and r = r{v). In all plots, the vertical 
axis is height, with the lower boundary at the piston's lowest point, and the upper 
boundary being the upper wall of the container. Solid line shows the position of the 
piston. The heavy dotted line shows the height j/cm of the center of mass. The upper 
and lower thin dotted lines show 2/cm ~t~ o'cm and Ucm — cFcm respcctivcly, where (Tcm 
is the standard deviation of the particle heights at a given instant. The majority of 
the particles are thus contained between the two thin dotted lines. 

vanishes because no particles hit the top of the container. The cluster bounces 
on the vibrating plate. As the piston velocity increases, the cluster flies higher 
and higher until it strikes the top. The pressure reaches a first maximum near 
V = 0.82 m/s (Fig. 11) where there is a first resonance between the time 
of flight of the cluster and the vibration period. These two times are such 
that the cluster lands on the piston just when the piston is at its maximum 
upwards velocity (Fig. 12b). Note that a large number of particles descend 
below the piston's maximum height. Then as the piston rises, it sweeps up all 
the particles, increasing the density of the cluster. On the other hand, as the 
cluster falls from the stationary wall, it slowly expands. 
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As the piston velocity increases further, the pressure decreases, because 
the cluster lands on the piston when it is near its maximum height (Fig. 12c). 
Then the pressure rises to a new maximum near V = AAf = 2.2 m/s (Fig. 11) 
because there is a second resonance: the cluster's flight time is equal to two 
piston cycles (Fig. 12d). Note that the frequencies of these two resonances 
(/ = 8.2 Hz and / = 21.6 Hz) are not related by a simple factor of 2. The 
reason for this is that the time of flight of the cluster depends in a complicated 
way on the piston velocity. Comparing Figs. 12b and 12d, one can sec that the 
velocity of the cluster when it leaves the piston is close to V, but the time it 
takes the cluster to fall back from the upper wall onto the vibrating is roughly 
independent of V . Note also that the collisions between the cluster and the 
piston or wall takes a finite time, as can be seen by comparing the times of 
the minima of ycm + CTcm, Vcm, and j/cm — in Fig. 12b. All of these efli'ects 
exclude a simple resonance condition for the pressure maxima. The resonant 
frequencies depend on g and h — ho. The resonant frequencies are proportional 
to g^/^ but decrease for as h — ho is increased. The resonant frequencies are 
not sensitive to changing n, as long as n is sufficiently large. 

At higher driving frequencies (Fig. 12e and Fig. 12f), the center of mass 
of the particles changes very little over one cycle, and there are no more 
resonances. In fact, setting g = changes Fig. 12f very little. The reason 
for this is that average particle kinetic energy is much larger than gh, so 
gravity is irrelevant. Higher order resonances, (where the cluster's flight time 
equals three or more vibration periods) are not observed. These resonances 
are probably not observed because there is not enough time for the particles 
to fall below the piston's highest point. Thus, the piston does not gather up 
all the particles into a single dense cluster, but rather acts to disperse the 
cluster. For example, imagine that the driving frequency is suddenly doubled 
in Fig. 12d just at the moment when the cluster starts its descent from the 
wall. When the cluster arrives at the piston, perhaps a half of the particles 
will be swept up by one vibration cycle, and a half in the next cycle. This is 
how energetic forcing breaks up clusters. 

5 Conclusions 

In Sec. 3, we brought simulations as close as possible to the experiments. To 
bring simulations into agreement with the experiments, it is essential to use 
a velocity dependent coefficient of restitution. It is especially important to 
take into account plastic deformations that cause the restitution coefficient 
to decrease rapidly with increasing impact velocity. Indeed, the restitution 
coefficient for strongly vibrated steel spheres is very far from the constant 
values of r = 0.95 or r = 0.9 that are often cited in simulations as typical for 
steel spheres. We also noted that it is very important to take into account the 
number of particle layers n. The dependence of the pressure P on the piston 
velocity V changes with n. It is not accurate to speak of "a" scaling exponent 
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for the pressure in terms of V: this exponent depends continuously on n, and 
docs not exist at high density (n > 3) under gravity, due to the clustering 
instabihty. 

In Sec. 4, we investigated this question of clustering. When there is no 
gravity, a loose cluster forms against the stationary wall opposite the vibrat- 
ing piston. Due to the velocity dependent restitution coefficient, this cluster 
acts as a buffer for fastest particles, and leads to a simple non-extensive scal- 
ing P cx N'^V. Unlike the case with gravity, the scaling holds for a wide 
range of parameters. We discussed the origin the scaling, and remarked that 
it should be observable in microgravity experiments. We also studied the effect 
of clustering in the presence of gravity. In this case, the clusters give a very 
complicated behavior, because there are interactions between the period of 
vibration and the flight time of the cluster. Finally, if the vibration frequency 
is strong enough, the cluster can be broken up, because the cluster-piston 
collision time becomes longer than vibration period. 
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